i=1∑nj=1∑nijgcd(i,j)
k=1∑ni=1∑nj=1∑nij[gcd(i,j)=k]
k=1∑nk3i=1∑⌊kn⌋j=1∑ij[gcd(i,j)=1]
k=1∑nk3d=1∑⌊kn⌋μ(d) d2i=1∑⌊kdn⌋j=1∑⌊kdn⌋ij
令sum2(n)=∑i=1ni2,sum3(n)=∑i=1ni3
利用小学奥数知识有:
sum2(n)=6n(n+1)(2n+1) , sum3(n)=4n2(n+1)2
上式为:
k=1∑nk3d=1∑⌊kn⌋μ(d) d2sum3(kdn)
记 T=kd , 则原式可继续化简为
T=1∑nk∣T∑k3μ(kT)(kT)2sum3(Tn)
T=1∑nsum3(Tn) T2k∣T∑kμ(kT)
后面是一个卷积形式: id∗μ=φ
T=1∑nsum3(Tn) T2φ(T)
现在只需求得 φ×id×id 的前缀和就可以整除分块了。
f(n)=φ(n)×n×n
g(n)=n×n
h(n)=d∣n∑f(d)∗g(dn)
=d∣n∑φ(d)×d×d∗dn∗dn
=n2d∣n∑φ(d)=n3
g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋)
S(n)=sum3(n)−d=2∑ng(d)S(⌊dn⌋)
这道题就完了。
#include <map>
#include <cstdio>
using namespace std;
const int MAXN = 8000000;
long long n;
int p , inv4 , inv6 , k , prime[ MAXN + 5 ] , phi[ MAXN + 5 ] , f[ MAXN + 5 ];
bool vis[ MAXN + 5 ];
void sieve( ) {
f[ 1 ] = phi[ 1 ] = 1;
for( int i = 2 ; i <= MAXN ; i ++ ) {
if( !vis[ i ] ) {
prime[ ++ k ] = i;
phi[ i ] = i - 1;
}
for( int j = 1 ; j <= k && 1ll * i * prime[ j ] <= MAXN ; j ++ ) {
vis[ i * prime[ j ] ] = 1;
if( i % prime[ j ] == 0 ) {
phi[ i * prime[ j ] ] = phi[ i ] * prime[ j ];
break;
}
else
phi[ i * prime[ j ] ] = phi[ i ] * ( prime[ j ] - 1 );
}
f[ i ] = ( f[ i - 1 ] + 1ll * i * i % p * phi[ i ] % p ) % p;
}
}
map< long long , int > Map;
int sum2( long long x ) {
x %= p;
return 1ll * x * ( x + 1 ) % p * ( 2 * x + 1 ) % p * inv6 % p;
}
int sum3( long long x ) {
x %= p;
return 1ll * x * x % p * ( x + 1 ) % p * ( x + 1 ) % p * inv4 % p;
}
int Sumphii( long long n ) {
if( n <= MAXN ) return f[ n ];
if( Map[ n ] ) return Map[ n ];
int Ans = sum3( n );
for( long long l = 2 , r ; l <= n ; l = r + 1 ) {
r = n / ( n / l );
Ans = ( Ans - 1ll * ( sum2( r ) - sum2( l - 1 ) ) % p * Sumphii( n / l ) % p ) % p;
}
return Map[ n ] = ( Ans + p ) % p;
}
int solve( long long n ) {
int Ans = 0;
for( long long l = 1 , r ; l <= n ; l = r + 1 ) {
r = n / ( n / l );
Ans = ( Ans + 1ll * ( Sumphii( r ) - Sumphii( l - 1 ) ) * sum3( n / l ) % p ) % p;
}
return ( Ans + p ) % p;
}
int Quick_pow( int x , int po ) {
int Ans = 1;
while( po ) {
if( po & 1 ) Ans = 1ll * Ans * x % p;
x = 1ll * x * x % p;
po /= 2;
}
return Ans;
}
int Inv( int x ) {
return Quick_pow( x , p - 2 );
}
signed main( ) {
scanf("%d %lld",&p,&n);
sieve( );
inv4 = Inv( 4 ) , inv6 = Inv( 6 );
printf("%d", solve( n ) );
return 0;
}